Estimating Concentration Response Function and Change-Point using Time-Course and Calibration Data

In this paper the problem of determining the functional relationship between time and the concentration of a chemical substance is studied. An intervention drug is administered on the experimental unit from which the chemical substance (specimen) is measured. This drug is hypothesized to cause a change in the concentration level of the chemical substance a certain lag-time after the intervention. However, the concentration value could not be directly measured, but rather a surrogate response can be measured. In the time-course study, this surrogate response is measured using different electrodes which possess varied behaviors. To utilize these surrogate measurements arising from the different electrodes (sensors), a calibration study is undertaken which measures the surrogate response for the different electrodes at known concentration levels. Based on the time-course and calibration data sets, a statistical procedure to estimate the signal function and the lag-time is proposed. Simulation studies indicate that the proposed procedure is able to reasonably recover the signal function and the lag-time. The procedure is then applied to the real data sets obtained during an analytical chemistry experiment.


Introduction
The general problem tackled in this paper is determining the time-behavior of the concentration of a chemical obtained from an experimental unit (e.u.) subjected to an intervention. Specifically, of interest is to study if a specific drug intervention causes a change in concentration. Concentration levels, however, are difficult to measure directly, especially when the e.u. could not be sacrificed, such as during a time-course study. Thus, the concentration levels are to be inferred indirectly through electric charge measurements which are stochastically related to the concentration levels. The measurement of the charge This work is licensed under Creative Commons Attribution 4.0 Licens is reliant on the use of an electrode. But different electrodes possess different behaviors when measuring the charge. Thus, to be able to use different electrodes to infer the concentration levels, a calibration study using these electrodes is performed to determine the relationship between concentration and charge for each of the electrodes.
The specific study that motivated the problem considered in this paper was performed in the laboratory of one of the authors (Dr. P. Hashemi) at the Department of Chemistry and Biochemistry, University of South Carolina. The chemical concentration of interest is that of serotonin, a substance that that has been implicated in affective disorders such as depression, and the e.u.'s are mice (in vivo). The intervention performed on the mice is the administration of a pharmacological agent or drug. The type of data sets obtained from the study are shown in Figures 1 & 2. Figure 1 shows charge measurements over time for five different electrodes, with panels 1 and 2 showing results for the two intervention agents or drugs: pargyline and GBR 12909, respectively. Pargyline inhibits serotonin metabolism while GBR 12909 is a dopamine reuptake inhibitor. These drugs or agents inhibit serotonin metabolism and dopamine reuptake. Figure 2 presents the calibration data for ten electrodes, together with fitted mean response curves. The five electrodes used in the time-course study for each of the two drugs came from these 10 electrodes utilized in the calibration study. For instance, the electrodes used in the time-course experiment with pargyline were the electrodes labeled 2, 3, 4, 5, and 10; whereas, for GBR 12909 the electrodes were those labeled 1, 6, 7, 8, and 9.
The goal in this study is to estimate the mean concentration response function over time based on the time course study and the calibration study. Another aspect in the time-course study is that the drug intervention was administered after an initial no-drug period. Thus, another important goal is to determine the lag-time post administering the drug after which the drug takes effect, if indeed the drug has an effect. Discipline-specific details of this motivating and focused application has been previously published [1]. We point out the statistical methodology developed in this paper and the ideas contained herein have the potential of being useful in other situations where calibrated measurements are obtained [2].

Mathematical setting
In this section we describe the postulated models for both the time-course study and the calibration study. For this purpose, let us suppose that we have a specimen from the e.u. We denote by X t the concentration level of the chemical at time t for this specimen. We denote by T the time at which the intervention is performed, that is, the time of administering the drug intervention. We shall postulate the stochastic model where ∆ is a lag-time and g(·) is a continuous function with g(t) = 0 for t ≤ 0. If smoothness is desired, we may impose the condition that g(·) is differentiable for t > 0. The time of intervention T is known, while the lag-time ∆, together with the regression coefficients γ 0 and γ 1 , and the function g(·)are unknown. The error is assumed to be white noise, with ε t having a normal distribution with mean zero and variance τ 2 , with τ also unknown.
As discussed in section 1, the X t could not be directly measured. If it could be directly measured and we possess the values of X t ′s, then we could estimate the model parameters, and also infer the change point T + ∆, or the lag-time ∆.
To enable the determination of the chemical concentration levels X t ′s, there are K possible electrodes that could be used for measuring the charge. At time t, n t charge measurements using different electrodes will be taken. Thus, denoting by Y the charge and by E the electrode type, at time t, n t pairs (Y tj , E tj ), j = 1, 2, …, n t , are obtained. The measured charge is affected by the chemical concentration level and the electrode type. The relationship between the charge and the concentration and electrode type is given by where β 0 , β 1 ,ξ 2 ,…ξ k ,η 2 ,…η k are unknown, and ∈ has a normal distribution with mean zero and variance σ 2 , with σ unknown. The function I(·) is the indicator function, taking a 1(0)-value depending on whether the argument is true(false). This is a linear model that incorporates an interaction between the concentration and the electrode type [3,4].
In order to estimate the model parameters in model (2), a calibration study is performed. See [5] for some review of calibration methods. In this study, we followed the classic inverse regression approach. First, we regress the response variable, Y, on predictor X, and estimate linear regression coefficients using least-squares method; thereafter, the value of an unknown X, is to be estimated given an observation of Y, by subtracting the estimated intercept and dividing by the estimated slope [6][7][8][9]. Inference for the calibration parameters, is not trivial because of the presence of a normally distributed estimated slope in the denominator, which causes the inverse estimator to have infinite variance [10]. In this paper, we use Delta method to construct approximate confidence intervals for the calibration parameters. A more conservative confidence interval approach based on inverting simultaneous tolerance intervals was proposed by Scheffe [9] in literature. An alternative approach to the problem is referred to as reverse regression, when X's are treated as the response and formally regressed on Y's (even though the X's are measured with negligible error). Krutchkoff [4] compared inverse and reverse regression using Monte Carlo simulations. Properties and limitations of the reverse estimators were studied by Williams [10] and Halperin [3].
In this calibration study, known levels of concentration are used, and the different electrodes are used to measure the charge. We denote by L the number of concentration levels, and these levels will be denoted by x 01 < x 02 < … < x 0L . At concentration level x 0l , all K electrode types are utilized, and for each electrode type, there are m charge measurements. Thus, for the l th concentration level, there are Km observations, and for the L concentration levels, there are N = LKm charge measurements. The data could be summarized as in Table   1 and pictorially depicted as in Figure 2.
With ξ 1 = η 1 = 0, the linear models governing these charge measurements are given by Y ijl = β 0 + ξ j + β 1 + n j x 0i + ∈ ijl (3) for i = 1,…, L; j = 1,…, K; l = 1,…, m, and with the ∈ ijl ′ s being independent and identically distributed (IID) and having a normal distribution with mean zero and variance σ 2 . These linear models could be fitted using object functions in a variety of statistical packages, such as the function lm or glm in the R statistical package [8].

Estimating parameters
Instead of using the calibration data representation presented in Table 1, for purposes of describing more concisely the estimators of parameters, we denote by Y = (Y 1 ,Y 2 ,…,Y N ) T the N × 1 vector of charge values. The design matrix is W, which is an N × 2K matrix whose ith row is , v Ki x i with x i the concentration level, and v ji = I E i = j , j = 2,…K, indicates whether the electrode type is j. With ∈ = ∈ 1 , ∈ 2 , …, ∈ N T denoting the error vector, the linear model could be written via Y = W θ + ∈ (4) where the 2K × 1 regression coefficient vector θ is θ = β 0 , β 1 , η 2 , …, η K , ξ 2 , …, ξ K T . (5) In this model, ∈ N N 0, σ 2 I N where I N is the N × N identity matrix.
The least-squares (LS) estimator of θ, which is also the maximum likelihood (ML) estimator, is given by (see, for instance, any linear theory book such as [6]) This is an unbiased estimator of θ. The error variance σ 2 is unbiasedly estimated by where H = W W T W −1 W T . An unbiased estimator of the covariance matrix of θ is Qiang et al. Page 4 The validation of the assumptions underlying this linear model could be performed graphically or via the global method in [7].

Calibrated concentration estimators
The calibration problem will now be described. Suppose that, at a given time, we are given a charge measurement for a specific electrode: say, (y 0 , e 0 ), where e 0 is the electrode type and y 0 is the charge on that specific electrode. More generally, suppose that M pairs of charge and electrode type are available at a given time: (y 01 , e 01 ), (y 02 , e 02 ),…, (y 0M , e 0M ).
Based on these observations, what is an estimate of the concentration level, and what is an approximate 100(1−α)% confidence interval for the concentration level?
By the Delta-Method, an estimate of the variance of x 0 y 0 , e 0 is V ar x y 0 , e 0 ; θ ≡ V 0 2 y 0 , e 0 = σ 2 b e0 y 0 , θ T W T W −1 b e0 y 0 , θ + 1 An approximate 100(1−α)% confidence interval for the concentration level, having observed a charge value of y 0 using electrode type e 0 , is given by Calibration based on many pairs (y 0 , e 0 ) Next we consider the situation where several charge measurements are taken using possibly different electrode types. Let y 0 , e 0 = y 0m , e 0m : m = 1, 2, …, M .
Here y 0m is the mth charge measurement which is obtained using electrode type e 0m . From the preceding subsection, based on this particular measurement we obtain an estimate of the concentration level, given by x 0 y 0 , e 0 ; θ which has an approximate estimated variance of V 0m 2 y 0m , e 0m whose expression is obtained via (11). The M estimates of the concentration level obtained for each of the elements in (y 0 , e 0 ) will not be independent of each other, owing to the fact that they will all depend on θ . In fact, through the Delta-Method, we could obtain their approximate covariance matrix.
However, for the sake of simplicity and practicality, we ignore the dependence among these M estimates. Under this simplified assumption, whose appropriateness will be examined later through simulation studies, we could then combine the M estimates by simply taking into account the possibly varying estimates of their variances. We recall the following well-known distribution theory result, which is easily proved using a Lagrange multiplier minimization approach.

Theorem 1
Let W 1 ,W 2 ,…,W m be independent random variables with common mean μ and respective variances τ 1 2 , τ 2 2 , …, τ m 2 . Among all linear combinations ∑ l = 1 m c l W l with ∑ l = 1 m c l = 1, so that the mean of the linear combination is still μ, the one with the smallest variance coincides with the choice of coefficients given by The variance of this optimal linear combination is Using this result, we obtain our combined estimate of the concentration level from the M estimates via x 0 y 0 , e 0 ; θ = ∑ m = 1 m c m x 0 y 0m , e 0m ; θ , where the weights are given by c m = 1/V 0m 2 y 0m , e 0m ∑ q = 1 M 1/V 0q 2 y 0q , e 0q , m = 1, 2, …M .

Author Manuscript
Author Manuscript Author Manuscript Author Manuscript

Time-course study
As mentioned in section, 1 the main goal of the study is to determine the behavior of the concentration of the chemical over a period of time, where at some point T during the monitoring period, a drug intervention is performed. We suppose that over the period [0,T*] charge measurements using possibly different electrode types are performed at specified times 0 ≤ t 1 < t 2 , < t 3 <…< t L ≤T*. The charge and electrode observations at time t 1 are given by (y 0l , e 0l ). Based on these observation vectors, we find an estimate of the concentration level at time t 1 to be x 0l y 0l , e 0l , together with its estimate of its variance V 0l 2 y 0l , e 0l ; θ . For each time t l , we could also construct the approximate 100(1− α)% confidence interval for the concentration value given by x 0l y 0l , e 0l ; θ ± t N − 2K; α/2 V 0l y 0l , e 0l .
Note, however, that these confidence intervals are not adjusted for multiplicity. On the other hand, based on the time course data and the resulting concentration estimates at each of the times of observations, given by t l , x 0l y 0l , e 0l ; θ , V 0l y 0l , e 0l : l = 1, 2, …, L , we could fit an appropriate regression curve that takes into account the possibly differing variability of the x 0l y 0l , e 0l ; θ ′ s. One of the simplest models that could be fitted to this data, under the hypothesis that upon the drug intervention at time T the concentration curve should change after a lag-time of ∆, specifies that x 0l = ω 0 + ω 1 t l + k 1 U l + k 2 U l 2 + V 0lεl , and U l = max 0, t 1 − T + Δ , l = 1, 2, …, L, and ε l has mean zero and variance τ 2 . Note that the τ 2 in this model is not the same as the τ 2 in the time-course model in (1). Here, the time of intervention T is known, whereas the lag-time ∆ is not known. For a specified ∆, this model is easily fitted using the lm command in R with the weights option enabled, which performs a weighted linear regression fit [8].
The coefficient of determination could be obtained and denoted by R 2 (∆). The coefficient of determination could be plotted with respect to possible values of ∆. The value of ∆ that maximizes this coefficient of determination R 2 (∆) is a plausible estimate for the lag-time ∆.
Thus, a possible estimator of the lag-time is This could be computed by fitting the model over a sequence of values of ∆ and searching for the value of ∆ that maximizes R 2 (∆). This computational method is the approach we implemented in this paper. The confidence interval described above at a given time point is appropriate if we only have data at one time point. However, in the time-course study, we actually obtain a series of concentration estimates for each of the time points considered. These pairs of time and concentration estimates are then used to estimate the time course model parameters. In the process of fitting the linear-parabolic curve, we could then also obtain pointwise confidence intervals for the concentration value at each of the time points, where we utilize the estimate of the error standard deviation. We surmise that the resulting point-wise confidence intervals is a better indicator of where the concentration values are.
In the real-data application of the procedure in section 6, this point-wise confidence interval approach will be presented.

A real-data illustration
This section provides a more detailed description of the statistical methods performed in the data analysis in some portions of paper [1]. The data sets were obtained in the Hashemi laboratory at University of South Carolina. We limit our consideration to the intervention drug pargyline, which was hypothesized to have an impact on the concentration level of serotonin, in contrast to the intervention drug GBR 12909. Figures 1-2 present the timecourse data and the fitted linear models. These figures are adapted with permission from [1] and are copyrighted from American Chemical Society. In the next section, we then present simulation studies to provide us some ideas on the properties of this procedure.
A time-course study was performed leading to the data set with charge measurements over time (from 0 minutes to 121 minutes), with the intervention drug pargyline (10 mg/kg) administered at 60 minutes, via interperitoneal injection. The time-course data with pargyline as the intervention drug is depicted in the first panel of Figure 1. The different symbols correspond to the 5 different electrodes (out of 10 possible electrodes) used in measuring the charges. A calibration study was also performed where, for known concentration levels, charge measurements were obtained for each of 10 possible electrodes. The charges of each electrode are measured under three concentration values (25, 50, and 100 nM) each with 4 replicates. The data set obtained from this calibration study are the solid points in Figure 2.
A linear model, as described in (2), was fitted using the calibration data. The fitted linear models, using the calibration data, for each electrode type are depicted as the lines in Figure  2 for each of the 10 electrodes. A summary of the estimated model parameters are in Table 2. These linear models with interaction terms provide excellent fits to the observed calibration data with coefficient of determinations all above 99%. The interaction effects are all significant.
Using these fitted linear models, given the charge measurements at each time point from the time course study, estimates of the concentration levels at each time point were obtained. This procedure yielded pairs of values of time and concentration estimates which are the solid points in Figure 3. A functional continuous model, as described in Section 5, was fitted to these pairs of time and concentration values using weighted regression via the lm command in the R environment. The estimate of the time-lag ∆ after which the drug takes effect is Δ = 3.32minutes. As mentioned earlier, this estimate was obtained by maximizing the coefficient of determination with respect to the possible values of ∆ [see discussion prior to (18)]. The resulting fitted linear-parabolic model whose equation is given in (18) is shown in Figure 3.  Table 3. Observe that the coefficient for the time effect is not significantly different from (p = .429) which is to be expected since without drug intervention the serotonin concentration is not expected to change. We also mention that the final fitted model has an R 2 equal to 95.59%, indicating an excellent fit of the linear-parabolic model for relating concentration to time for this study.

Simulation studies
In order to study the properties of the procedure described above for estimating the lag-time parameter ∆, we performed a computer simulation study. All simulation programs were coded in R and function objects in R such as lm, rnorm, etc. were utilized. The main purpose of this study was to determine if we are able to effectively estimate the lag-time parameter ∆ using a time-course data and a calibration data. In the simulation, we generate the time-course data using the model Where Z tm 2 ′ s are IID random variables. The parameter values were set to β 0 = 1, β 1 = 1.5, ξ = (0,1,2), η = (0,.5,1.5). Thus, this leads to the 'observed time-course data' given by t, y tm , m = 1, 2, 3; t = 0, 1, 2, …120 .
A calibration experiment was also simulated. In this calibration experiment, five specific concentration values the goal is to see if we could reasonably recover the signal function given by g t = E X t = γ 0 + γ 1 t + k 1 max 0, t − T + Δ + k 2 max 0, t − T + Δ 2 (22) and obtain good estimates of the model parameters using the procedure described in the preceding sections, in particular to obtain a reasonable estimate of ∆.
We performed the above-described simulation experiment using 5000 replications. The estimates of ∆ are provided in the frequency histogram in Figure 4. In searching for the optimizing ∆, we used increments of 0:1, so that the possible values of Δ contains only one decimal digit. The mean of these Δ ′ s is 5:03828 with a standard error of 0:4990. Thus, this mean is quite close to the true value of ∆ which is 5.0. Regarding the other model parameters pertaining to the time-course portion, we summarize the results in Table 4. Histograms of these 5000 estimates for the four model parameters are depicted in Figure 5. Observe that the means of these estimates for each of the regression parameters are close to their true values.
Another parameter in the time-course model is τ, which is the standard deviation of the error component. When we fit the linear-parabolic model on the time-course and calibration data, an estimate of the error standard deviation is also obtained. However, this is not an estimate of τ since this estimates a larger value than τ owing to the additional error contamination arising from the calibration study and the charge measurements at each time and for each electrode type. For instance, in the simulation study, the histogram of the 5000 estimates of the standard deviation of the error term provided in Figure 6 all exceed the true value of τ = 10. The mean of these standard deviation estimates is 23.56979 with a standard deviation of 2.723542. On the basis of this modest simulation study, the statistical procedure developed in earlier sections using the calibration and time-course data sets appears successful in recovering the functional relationship between time and concentration and is also able to infer properly about the lag-time at which the drug starts to become effective.
In order to further study the properties of the statistical procedure on the real data, we performed another simulation study using the estimated model parameters based on the data from the Hashemi laboratory. The time course data was generated using model (20) where the parameter values were set identical to the estimated values in Table 2, and with x 1 = 25, x 2 = 50, x = 100, and Z jmk 3 ′ s being IID from a N(0, 0.095 2 ) distribution. Using the time-course data and the calibration data, it is of interest to see if we could reasonably recover the signal function in (22) and also obtain good estimates of the model parameters.
The resulting 5000 estimates of ∆ are provided in the frequency histogram in Figure 7.
The mean of these 5000 Δ ′ s is 3:3008 with a standard error of 1:9341. The mean is indeed very close to the true value which is ∆ = 3.3, indicating the possible unbiasedness of the estimator, though of course a simulation study could not establish this theoretical property. Note, however, that there were some negative estimates that arose and this could be a consequence of a wrong (we hypothesize larger values) specification of the error variances in the models. Details of the estimates for the other parameters are summarized in Table 5.
Histograms of the estimates for these four model parameters are depicted in Figure 8.

Concluding Remarks
Motivated by a study in an analytical chemistry laboratory dealing with the concentration of the chemical substance serotonin, we developed in this paper a statistical procedure for estimating the functional relationship between time and concentration level of a substance, together with the change point, based on a time-course study that measures a concentrationsurrogate variable (the charge in this application) and data from a calibration study. The novel aspect of this procedure is the presence of several measuring electrodes which have their unique behaviors when utilized to measure the surrogate variable (the charge in application). Simulation studies were performed to examine the properties of the proposed procedures, and the simulation results indicate that the procedure is able to reasonably recover the signal function and also the change-point in the signal function.
The procedure is also applied to the time-course data and the calibration data from the focus application obtained in the analytical chemistry experiment. Further theoretical studies and simulations will be desired to examine more carefully the procedure especially when applied to more complex settings such as when the g (·) portion in the signal function in equation (1) is not parabolic or when it is non-parametrically specified. Of interest for further research are situations where the error distributions in the regression models are not normally-distributed, which may entail the use of non-parametric regression estimation methods.    Frequency histogram of the estimates of from the simulation study. The true value was and there were 5000 replications.  Frequency histogram of the estimates of from the simulation study with 5000 replications. The true value was ∆ = 3.3.  Format of the calibration data from the calibration study.

Electrode Type
Concentration Levels     Means and standard deviations of the 5000 estimates obtained in the simulation study of the four model parameters  Means and standard deviations of the 5000 estimates obtained in the simulation study of the four model parameters